C
C  Returns the day-of-the-week (Sunday=1) of the given daten
C
      integer*4 function dayofwk (y, m, d)
C
      implicit none
      integer*4 y, m, d, tmp1, tmp2
      automatic tmp1, tmp2

      tmp1 = m + 10
      tmp2 = y + (m - 14) / 12
      dayofwk = mod ((13 * (tmp1 - tmp1 / 13 * 12) - 1) / 5 +
     +          d + 77 + 5 * (tmp2 - tmp2 / 100 * 100) / 4 +
     +          tmp2 / 400 - tmp2 / 100 * 2, 7) + 1
      return
      end
C
C  Return the day-of-the-year corresponding to the given date
C
      integer*4 function dayofyr (y, m, d)
C
      implicit none
      integer*4 y, m, d, itmp
      automatic itmp

      itmp = (m + 10) / 13
      dayofyr = 3055 * (m + 2) / 100 - itmp * 2 - 91 +
     +          (1 - (y - y / 4 * 4 + 3) / 4 +
     +          (y - y / 100 * 100 + 99) / 100 -
     +          (y - y / 400 * 400 + 399) / 400) *
     +          itmp + d
      return
      end
C
C  Compute month and day given year and day-of-the-year
C
      subroutine calndr1 (dayofyr, y, m, d)
C
      implicit none
      integer*4 dayofyr, y, m, d, t1, t2, t3, t4
      automatic t1, t2, t3, t4

      t1 = 0
      if (mod (y, 4) .eq. 0) t1 = 1
      t2 = 0
      if ((mod (y, 400) .eq. 0) .or.
     +    (mod (y, 100) .ne. 0)) t2 = t1
      t1 = 0
      if (dayofyr .gt. t2 + 59) t1 = 2 - t2
      t3 = dayofyr + t1
      t4 = ((t3 + 91) * 100) / 3055
      d = (t3 + 91) - (t4 * 3055) / 100
      m = t4 - 2
      return
      end
C
C  Return the day-of-the-year of the first day of a given year and month
C
      function first_mon (y, m)
C
      implicit none
      integer*4 first_mon, y, m, yy, dayofyr, dayofwk
      automatic yy

      external dayofyr, dayofwk
      yy = y
      first_mon = dayofyr (yy, m, 1) + mod (9 - dayofwk (yy, m, 1), 7)
      return
      end
C
C  Compute month <m> and day-of-the-month <d> for Easter Sunday
C  in a given year (yr>.
C
      subroutine easter (yr, m, d)
      implicit none
      integer*4 yr, t1, t2, t3, t4, t5, m, d
      automatic t1, t2, t3, t4, t5

      t1 = mod(yr,19)
      t2 = yr/100
      t3 = mod(yr,100)
      t4 = mod(19*t1+t2-t2/4-((t2-(t2+8)/25+1)/3)+15,30)
      t5 = mod(32+2*mod(t2,4)+2*(t3/4)-t4-mod(t3,4),7)
      t2 = t4+t5-7*((t1+11*t4+22*t5)/451)+114
      m = t2/31
      d = mod(t2,31)+1
      return
      end
C
C   Returns Julian day number for any Gregorian Calendar date *)
C   GIVEN:   <y> - year <m> - month and <d> - day
C
      integer*4 function julday (y, m, d)
      implicit none
      integer*4 y, m, d, t1, t2, t3, t4
      automatic t1, t2, t3, t4
C
      if (m .gt. 2) then
         t1 = m - 3
         t2 = y
      else
         t1 = m + 9
         t2 = y - 1
      endif
      t3 = t2 / 100
      t4 = mod (t2, 100)
      julday = (146097 * t3) / 4 + (1461 * t4) / 4 +
     +         (153 * t1 + 2) / 5 + d + 1721119
      return
      end

C
C  Computes calendar date from a given Julian day number <julday> for
C  any valid Gregorian calendar date <y>, <m>, <d>
C
      subroutine calndrj (julday, y, m, d)
      implicit none
      integer*4 julday, y, m, d, t1, t2
      automatic t1, t2

      t1 = julday - 1721119 
      y = (4 * t1 - 1) / 146097 
      t1 = mod (4 * t1 - 1, 146097)
      t2 = t1 / 4 
      t1 = (4 * t2 + 3) / 1461 
      t2 = mod (4 * t2 + 3, 1461)
      t2 = (t2 + 4) / 4 
      m = (5 * t2 - 3) / 153 
      t2 = mod (5 * t2 - 3, 153)
      d = (t2 + 5) / 5 
      y = 100 * y + t1 
      if (m .lt. 10) then
        m = m + 3
      else
        m = m - 9 
        y = y + 1
      endif
      return
      end
